rm(list=ls())
setwd("F:/n2mosurgery")
library(table1)
library(survival)
library(survminer)
library(rms)
library(pec)
library(survivalROC)
library(plyr)
library(regplot)
library(riskRegression)
library(pROC)
library(dplyr)
library(MatchIt)
library(tidyr)
#install.packages("pec")
library(pec)
data<-read.csv("1.csv",header=T,sep=",")
seer=data

seer$Age<-factor(seer$Age,labels = c("<60","60-69","70-79",">=80"))
seer$Race<-factor(seer$Race,labels=c("Black","White","Other","Unknow"))
seer$Gender<-factor(seer$Gender,labels=c("Female","Male"))
seer$Location<-factor(seer$Location,labels=c("RUL","RML","RLL","LUL","LLL","others"))
seer$Grade<-factor(seer$Grade,labels=c("I","II","III","IV","Unknow"))
seer$T<-factor(seer$T,labels=c("T1","T2","T3","T4"))
# seer$N<-factor(seer$N,labels=c("N0","N1","N2","N3"))
seer$Surgery<-factor(seer$Surgery,labels=c("No","Yes"))
seer$SRLNR<-factor(seer$SRLNR,labels =c("None","<4",">=4","Others"))
seer$Radiation<-factor(seer$Radiation,labels=c("No","Yes"))
seer$Chemotherapy<-factor(seer$Chemotherapy,labels=c("No","Yes"))
set.seed(888)
tr<-sample(nrow(seer),0.7*nrow(seer))
train<-seer[tr,]
test<-seer[-tr,]
seer$group<-0
seer$group[-tr]<-1
seer$group<-factor(seer$group,labels=c("train","test"))
fit<-survfit(Surv(Survival,Status)~group,data=seer)
summary(fit,times=c(12,36,60))

pvalue <- function(x, ...) {
  y <- unlist(x)
  g <- factor(rep(1:length(x), times=sapply(x, length)))
  if (is.numeric(y)) {
    p <- t.test(y ~ g)$p.value #数值型数据用t-test(两组比较)
  } else {
    p <- chisq.test(table(y, g))$p.value #因子型数据用卡方
  }
  c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}

table1(~Age+Race+Gender+Location+Grade+T+Surgery+SRLNR+Radiation+Chemotherapy|group,data=seer,extra.col=list("P-value"=pvalue),overall=F,topclass="Rtable1-zebra")

agecox<-coxph(Surv(Survival,Status)~Age,data=train)
agehr<-summary(agecox)
HRage<- round(agehr$coefficients[,2],2)
CI5age <-round(agehr$conf.int[,3],2)
CI95age <-round(agehr$conf.int[,4],2)
PValueage <- round(agehr$coefficients[,5],3)
mul_CIage<-paste(CI5age,'-',CI95age)
agehrdata<- data.frame("HR"=HRage,"CI"=mul_CIage, "P"=PValueage)
agehrdata$a<- "("; 
agehrdata$c<- ")"
agehrdata<- agehrdata[,c("HR","a","CI","c","P")]
agehrdata<-unite(agehrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)

racecox<-coxph(Surv(Survival,Status)~Race,data=train)
racehr<-summary(racecox)
HRrace<- round(racehr$coefficients[,2],2)
CI5race <-round(racehr$conf.int[,3],2)
CI95race <-round(racehr$conf.int[,4],2)
PValuerace <- round(racehr$coefficients[,5],3)
mul_CIrace<-paste(CI5race,'-',CI95race)
racehrdata<- data.frame("HR"=HRrace,"CI"=mul_CIrace, "P"=PValuerace)
racehrdata$a<- "("; 
racehrdata$c<- ")"
racehrdata<- racehrdata[,c("HR","a","CI","c","P")]
racehrdata<-unite(racehrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)


gendercox<-coxph(Surv(Survival,Status)~Gender,data=train)
genderhr<-summary(gendercox)
HRgender<- round(genderhr$coefficients[,2],2)
CI5gender <-round(genderhr$conf.int[,3],2)
CI95gender <-round(genderhr$conf.int[,4],2)
PValuegender <- round(genderhr$coefficients[,5],3)
mul_CIgender<-paste(CI5gender,'-',CI95gender)
genderhrdata<- data.frame("HR"=HRgender,"CI"=mul_CIgender, "P"=PValuegender)
genderhrdata$a<- "("; 
genderhrdata$c<- ")"
genderhrdata<- genderhrdata[,c("HR","a","CI","c","P")]
genderhrdata<-unite(genderhrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)

locationcox<-coxph(Surv(Survival,Status)~Location,data=train)
locationhr<-summary(locationcox)
HRlocation<- round(locationhr$coefficients[,2],2)
CI5location <-round(locationhr$conf.int[,3],2)
CI95location <-round(locationhr$conf.int[,4],2)
PValuelocation <- round(locationhr$coefficients[,5],3)
mul_CIlocation<-paste(CI5location,'-',CI95location)
locationhrdata<- data.frame("HR"=HRlocation,"CI"=mul_CIlocation, "P"=PValuelocation)
locationhrdata$a<- "("; 
locationhrdata$c<- ")"
locationhrdata<- locationhrdata[,c("HR","a","CI","c","P")]
locationhrdata<-unite(locationhrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)

gradecox<-coxph(Surv(Survival,Status)~Grade,data=train)
gradehr<-summary(gradecox)
HRgrade<- round(gradehr$coefficients[,2],2)
CI5grade <-round(gradehr$conf.int[,3],2)
CI95grade <-round(gradehr$conf.int[,4],2)
PValuegrade <- round(gradehr$coefficients[,5],3)
mul_CIgrade<-paste(CI5grade,'-',CI95grade)
gradehrdata<- data.frame("HR"=HRgrade,"CI"=mul_CIgrade, "P"=PValuegrade)
gradehrdata$a<- "("; 
gradehrdata$c<- ")"
gradehrdata<- gradehrdata[,c("HR","a","CI","c","P")]
gradehrdata<-unite(gradehrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)

tcox<-coxph(Surv(Survival,Status)~T,data=train)
thr<-summary(tcox)
HRt<- round(thr$coefficients[,2],2)
CI5t <-round(thr$conf.int[,3],2)
CI95t <-round(thr$conf.int[,4],2)
PValuet <- round(thr$coefficients[,5],3)
mul_CIt<-paste(CI5t,'-',CI95t)
thrdata<- data.frame("HR"=HRt,"CI"=mul_CIt, "P"=PValuet)
thrdata$a<- "("; 
thrdata$c<- ")"
thrdata<- thrdata[,c("HR","a","CI","c","P")]
thrdata<-unite(thrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)

surgerycox<-coxph(Surv(Survival,Status)~Surgery,data=train)
surgeryhr<-summary(surgerycox)
HRsurgery<- round(surgeryhr$coefficients[,2],2)
CI5surgery <-round(surgeryhr$conf.int[,3],2)
CI95surgery <-round(surgeryhr$conf.int[,4],2)
PValuesurgery <- round(surgeryhr$coefficients[,5],3)
mul_CIsurgery<-paste(CI5surgery,'-',CI95surgery)
surgeryhrdata<- data.frame("HR"=HRsurgery,"CI"=mul_CIsurgery, "P"=PValuesurgery)
surgeryhrdata$a<- "("; 
surgeryhrdata$c<- ")"
surgeryhrdata<- surgeryhrdata[,c("HR","a","CI","c","P")]
surgeryhrdata<-unite(surgeryhrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)


LNRcox<-coxph(Surv(Survival,Status)~SRLNR,data=train)
LNRhr<-summary(LNRcox)
HRLNR<- round(LNRhr$coefficients[,2],2)
CI5LNR <-round(LNRhr$conf.int[,3],2)
CI95LNR <-round(LNRhr$conf.int[,4],2)
PValueLNR <- round(LNRhr$coefficients[,5],3)
mul_CILNR<-paste(CI5LNR,'-',CI95LNR)
LNRhrdata<- data.frame("HR"=HRLNR,"CI"=mul_CILNR, "P"=PValueLNR)
LNRhrdata$a<- "("; 
LNRhrdata$c<- ")"
LNRhrdata<- LNRhrdata[,c("HR","a","CI","c","P")]
LNRhrdata<-unite(LNRhrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)

radiationcox<-coxph(Surv(Survival,Status)~Radiation,data=train)
radiationhr<-summary(radiationcox)
HRradiation<- round(radiationhr$coefficients[,2],2)
CI5radiation <-round(radiationhr$conf.int[,3],2)
CI95radiation <-round(radiationhr$conf.int[,4],2)
PValueradiation <- round(radiationhr$coefficients[,5],3)
mul_CIradiation<-paste(CI5radiation,'-',CI95radiation)
radiationhrdata<- data.frame("HR"=HRradiation,"CI"=mul_CIradiation, "P"=PValueradiation)
radiationhrdata$a<- "("; 
radiationhrdata$c<- ")"
radiationhrdata<- radiationhrdata[,c("HR","a","CI","c","P")]
radiationhrdata<-unite(radiationhrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)


chemocox<-coxph(Surv(Survival,Status)~Chemotherapy,data=train)
chemohr<-summary(chemocox)
HRchemo<- round(chemohr$coefficients[,2],2)
CI5chemo <-round(chemohr$conf.int[,3],2)
CI95chemo <-round(chemohr$conf.int[,4],2)
PValuechemo <- round(chemohr$coefficients[,5],3)
mul_CIchemo<-paste(CI5chemo,'-',CI95chemo)
chemohrdata<- data.frame("HR"=HRchemo,"CI"=mul_CIchemo, "P"=PValuechemo)
chemohrdata$a<- "("; 
chemohrdata$c<- ")"
chemohrdata<- chemohrdata[,c("HR","a","CI","c","P")]
chemohrdata<-unite(chemohrdata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)
hrdata<-rbind(agehrdata,racehrdata,genderhrdata,locationhrdata,gradehrdata,thrdata,surgeryhrdata,LNRhrdata,radiationhrdata,chemohrdata)
write.csv(hrdata,file="hrdata.csv")
print(genderhrdata)


traincox<-coxph(Surv(Survival,Status)~Age+Gender+Location+T+Surgery+SRLNR+Radiation+Chemotherapy,data=train)
mulHR<-summary(traincox) 
mul_HR<- round(mulHR$coefficients[,2],2) 
mul_PValue<- round(mulHR$coefficients[,5],4) 
mul_CI1<-round(mulHR$conf.int[,3],2)
mul_CI2<-round(mulHR$conf.int[,4],2)
mul_CI<-paste(mul_CI1,'-',mul_CI2)
muldata<- data.frame("HR"=mul_HR,"CI"=mul_CI, "P"=mul_PValue)
muldata$a<- "("; 
muldata$c<- ")"
muldata<- muldata[,c("HR","a","CI","c","P")]
muldata<-unite(muldata, "HR (95%CI)", c(HR, a, CI, c), sep = "", remove=T)
write.csv(muldata,file="mulhrdata.csv")

# y<- Surv(time=train$Survival,event=train$Status==1)
# 
# Uni_cox_model<- function(x){
#   FML <- as.formula(paste0 ("y~",x))
#   cox<- coxph(FML,data=train)
#   cox1<-summary(cox)
#   HR <- round(cox1$coefficients[,2],2)
#   PValue <- round(cox1$coefficients[,5],3)
#   CI5 <-round(cox1$conf.int[,3],2)
#   CI95 <-round(cox1$conf.int[,4],2)
#   Uni_cox_model<- data.frame('Characteristics' = x,
#                              'HR' = HR,
#                              'CI5' = CI5,
#                              'CI95' = CI95,
#                              'p' = PValue)
#   return(Uni_cox_model)}  
# names(train)
# variable.names<- colnames(train)[c(1:10)] 
# Uni_cox <- lapply(variable.names, Uni_cox_model)
# Uni_cox<- ldply(Uni_cox,data.frame)
# Uni_cox$CI<-paste(Uni_cox$CI5,'-',Uni_cox$CI95)
# Uni_cox<-Uni_cox[,-3:-4]
# 
# 
# Uni_cox$Characteristics[Uni_cox$p<0.05]
# 
# 
# mul_cox_model<- as.formula(paste0 ("y~",
#                                    paste0(Uni_cox$Characteristics[Uni_cox$p<0.05],collapse = "+")))
# mul_cox<-coxph(mul_cox_model,data=train)
# cox4<-summary(mul_cox) 
# mul_HR<- round(cox4$coefficients[,2],2) 
# mul_PValue<- round(cox4$coefficients[,5],4) 
# mul_CI1<-round(cox4$conf.int[,3],2)
# mul_CI2<-round(cox4$conf.int[,4],2)
# mul_CI<-paste(mul_CI1,'-',mul_CI2)
# mul_cox1<- data.frame("HR"=mul_HR,"CI"=mul_CI, "P"=mul_PValue)
# 
# library(survcomp)
# install.packages("survcomp")
# cox<-coxph(Surv(Survival,Status)~Age+Gender+Location+T+Surgery+SRLNR+Radiation+Chemotherapy,data=train)
# cindextrain<-concordance.index(pedict(fit),
#                                surv.time=train$Survival, surv.event=train$Status, method="noether")
#
summary(cox)
c_index<-survConcordance(Surv(train$Survival,train$Status)~predict(cox,train))$concordance
print(c_index)
c_index<-survConcordance(Surv(test$Survival,test$Status)~predict(cox,test))$concordance
# 
# print(c_index)
# c_index
# 计算总分
# dd<-datadist(train)
# option<-options(datadist="dd")
# coxm<-cph(Surv(Survival,Status)~Age+Gender+Location+T+Surgery+SRLNR+Radiation+Chemotherapy,data=train,surv = T,x=T,y=T)
# surv<-Survival(coxm)
# nom<-nomogram(coxm,fun=list(function(x)surv(12,x),function(x)surv(36,x),
#                             function(x)surv(60,x)),
#               lp=T,funlabel = c("12","36","60"))
# plot((nom),xfrac=.3)
# install.packages("nomogramFormula")
# library(nomogramFormula)
# options(option)
# results<-formula(nomogram=nom)
# points<-points_cal(formula=results$formula,lp=coxm$linear.predictors)
# print(points)
# summary(cox)


regplot(traincox,plots=c("boxplot","boxes"),observation=TRUE,failtime=c(12,36,60),title = "nomogram",points=TRUE,droplines=TRUE)


cox<-coxph(Surv(Survival,Status)~Age+Gender+Location+Grade+T+Surgery+SRLNR+Radiation+Chemotherapy,data=train,x=TRUE,y=TRUE)
dd<-datadist(train)
options(dataist="dd")
#train$lp<-predict(cox,type="risk",newdata =train)
train$lp1<-1-predictSurvProb(cox,newdata=train,times=12)
train$lp2<-1-predictSurvProb(cox,newdata=train,times=36)
train$lp3<-1-predictSurvProb(cox,newdata=train,times=60)
roctrain1y=survivalROC(Stime=train$Survival,
                  status=train$Status,
                  marker=train$lp1,
                  cut.values =1-0.627,
                  predict.time=12,method="KM")
roctrain3y=survivalROC(Stime=train$Survival,
                  status=train$Status,
                  marker=train$lp2,
                  cut.values = 1-0.326,
                  predict.time=36,method="KM")
roctrain5y=survivalROC(Stime=train$Survival,
                  status=train$Status,
                  marker=train$lp3,
                  cut.values = 1-0.205,
                  predict.time=60,method="KM")
plot(roctrain1y$FP,roctrain1y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctrain3y$FP,roctrain3y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctrain5y$FP,roctrain5y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.4,c(paste("12 mothons AUC:",round(roctrain1y$AUC,3)),
                 paste("36 mothons AUC:",round(roctrain3y$AUC,3)),
                 paste("60 mothons AUC:",round(roctrain5y$AUC,3))),
       x.intersp=0.1, y.intersp=0.5,
       lty= 1 ,lwd= 2,col=c("red","green","blue"),
       bty = "n",
       seg.len=0.5,cex=1)
roctrain5y$TP
1-roctrain5y$FP
Youde<-roctrain5y$TP-roctrain5y$FP
print(max(Youde))
print(which.max(Youde))
1-roctrain5y$cut.values[481]
roctrain1y
cox

#test$lp<-predict(cox,type="risk",newdata = test)
test$lp1<-1-predictSurvProb(cox,newdata=test,times=12)
test$lp2<-1-predictSurvProb(cox,newdata=test,times=36)
test$lp3<-1-predictSurvProb(cox,newdata=test,times=60)
roctest1y=survivalROC(Stime=test$Survival,
                status=test$Status,
                marker=test$lp1,
                cut.values =1-0.627,
                predict.time=12,method="KM")
roctest3y=survivalROC(Stime=test$Survival,
                  status=test$Status,
                  marker=test$lp2,
                  cut.values = 1-0.326,
                  predict.time=36,method="KM")
roctest5y=survivalROC(Stime=test$Survival,
                  status=test$Status,
                  marker=test$lp3,
                  cut.values = 1-0.205,
                  predict.time=60,method="KM")
plot(roctest1y$FP,roctest1y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctest3y$FP,roctest3y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctest5y$FP,roctest5y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.4,c(paste("12 mothons AUC:",round(roctest1y$AUC,3)),
                 paste("36 mothons AUC:",round(roctest3y$AUC,3)),
                 paste("60 mothons AUC:",round(roctest5y$AUC,3))),
       x.intersp=0.1, y.intersp=0.5,
       lty= 1 ,lwd= 2,col=c("red","green","blue"),
       bty = "n",
       seg.len=0.5,cex=1)
roctest5y$TP
1-roctest5y$FP


coxtrainage<-coxph(Surv(Survival,Status)~Age,data=train)
train$coxtrainagelp<-predict(coxtrainage,type="risk",newdata =train)
roctrainage1y=survivalROC(Stime=train$Survival,
                          status=train$Status,
                          marker=train$coxtrainagelp,
                          predict.time=12,method="KM")
roctrainage3y=survivalROC(Stime=train$Survival,
                          status=train$Status,
                          marker=train$coxtrainagelp,
                          predict.time=36,method="KM")
roctrainage5y=survivalROC(Stime=train$Survival,
                          status=train$Status,
                          marker=train$coxtrainagelp,
                          predict.time=60,method="KM")

coxtraingender<-coxph(Surv(Survival,Status)~Gender,data=train)
train$coxtraingenderlp<-predict(coxtraingender,type="risk",newdata =train)
roctraingender1y=survivalROC(Stime=train$Survival,
                             status=train$Status,
                             marker=train$coxtraingenderlp,
                             predict.time=12,method="KM")
roctraingender3y=survivalROC(Stime=train$Survival,
                             status=train$Status,
                             marker=train$coxtraingenderlp,
                             predict.time=36,method="KM")
roctraingender5y=survivalROC(Stime=train$Survival,
                             status=train$Status,
                             marker=train$coxtraingenderlp,
                             predict.time=60,method="KM")


coxtrainlocation<-coxph(Surv(Survival,Status)~Location,data=train)
train$coxtrainlocationlp<-predict(coxtrainlocation,type="risk",newdata =train)
roctrainlocation1y=survivalROC(Stime=train$Survival,
                               status=train$Status,
                               marker=train$coxtrainlocationlp,
                               predict.time=12,method="KM")
roctrainlocation3y=survivalROC(Stime=train$Survival,
                               status=train$Status,
                               marker=train$coxtrainlocationlp,
                               predict.time=36,method="KM")
roctrainlocation5y=survivalROC(Stime=train$Survival,
                               status=train$Status,
                               marker=train$coxtrainlocationlp,
                               predict.time=60,method="KM")

coxtraint<-coxph(Surv(Survival,Status)~T,data=train)
train$coxtraintlp<-predict(coxtraint,type="risk",newdata =train)
roctraint1y=survivalROC(Stime=train$Survival,
                        status=train$Status,
                        marker=train$coxtraintlp,
                        predict.time=12,method="KM")
roctraint3y=survivalROC(Stime=train$Survival,
                        status=train$Status,
                        marker=train$coxtraintlp,
                        predict.time=36,method="KM")
roctraint5y=survivalROC(Stime=train$Survival,
                        status=train$Status,
                        marker=train$coxtraintlp,
                        predict.time=60,method="KM")


coxtrainsurgery<-coxph(Surv(Survival,Status)~Surgery,data=train)
train$coxtrainsurgerylp<-predict(coxtrainsurgery,type="risk",newdata =train)
roctrainsurgery1y=survivalROC(Stime=train$Survival,
                              status=train$Status,
                              marker=train$coxtrainsurgerylp,
                              predict.time=12,method="KM")
roctrainsurgery3y=survivalROC(Stime=train$Survival,
                              status=train$Status,
                              marker=train$coxtrainsurgerylp,
                              predict.time=36,method="KM")
roctrainsurgery5y=survivalROC(Stime=train$Survival,
                              status=train$Status,
                              marker=train$coxtrainsurgerylp,
                              predict.time=60,method="KM")


coxtrainrlnr<-coxph(Surv(Survival,Status)~SRLNR,data=train)
train$coxtrainrlnrlp<-predict(coxtrainrlnr,type="risk",newdata =train)
roctrainrlnr1y=survivalROC(Stime=train$Survival,
                           status=train$Status,
                           marker=train$coxtrainrlnrlp,
                           predict.time=12,method="KM")
roctrainrlnr3y=survivalROC(Stime=train$Survival,
                           status=train$Status,
                           marker=train$coxtrainrlnrlp,
                           predict.time=36,method="KM")
roctrainrlnr5y=survivalROC(Stime=train$Survival,
                           status=train$Status,
                           marker=train$coxtrainrlnrlp,
                           predict.time=60,method="KM")


coxtrainradiation<-coxph(Surv(Survival,Status)~Radiation,data=train)
train$coxtrainradiationlp<-predict(coxtrainradiation,type="risk",newdata =train)
roctrainradiation1y=survivalROC(Stime=train$Survival,
                                status=train$Status,
                                marker=train$coxtrainradiationlp,
                                predict.time=12,method="KM")
roctrainradiation3y=survivalROC(Stime=train$Survival,
                                status=train$Status,
                                marker=train$coxtrainradiationlp,
                                predict.time=36,method="KM")
roctrainradiation5y=survivalROC(Stime=train$Survival,
                                status=train$Status,
                                marker=train$coxtrainradiationlp,
                                predict.time=60,method="KM")


coxtrainchemotherapy<-coxph(Surv(Survival,Status)~Chemotherapy,data=train)
train$coxtrainchemotherapylp<-predict(coxtrainchemotherapy,type="risk",newdata =train)
roctrainchemotherapy1y=survivalROC(Stime=train$Survival,
                                   status=train$Status,
                                   marker=train$coxtrainchemotherapylp,
                                   predict.time=12,method="KM")
roctrainchemotherapy3y=survivalROC(Stime=train$Survival,
                                   status=train$Status,
                                   marker=train$coxtrainchemotherapylp,
                                   predict.time=36,method="KM")
roctrainchemotherapy5y=survivalROC(Stime=train$Survival,
                                   status=train$Status,
                                   marker=train$coxtrainchemotherapylp,
                                   predict.time=60,method="KM")

# 训练集一年各变量比较图
plot(roctrain1y$FP,roctrain1y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctrainage1y$FP,roctrainage1y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctraingender1y$FP,roctraingender1y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
lines(roctrainlocation1y$FP,roctrainlocation1y$TP,type="l",col="yellow",xlim=c(0,1),ylim = c(0,1))
lines(roctraint1y$FP,roctraint1y$TP,type="l",col="orange",xlim=c(0,1),ylim = c(0,1))
lines(roctrainsurgery1y$FP,roctrainsurgery1y$TP,type="l",col="purple",xlim=c(0,1),ylim = c(0,1))
lines(roctrainrlnr1y$FP,roctrainrlnr1y$TP,type="l",col="gray",xlim=c(0,1),ylim = c(0,1))
lines(roctrainradiation1y$FP,roctrainradiation1y$TP,type="l",col="pink",xlim=c(0,1),ylim = c(0,1))
lines(roctrainchemotherapy1y$FP,roctrainchemotherapy1y$TP,type="l",col="skyblue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.6,c(paste("Nomogram AUC:",round(roctrain1y$AUC,3)),
                 paste("Age AUC:",round(roctrainage1y$AUC,3)),
                 paste("Gender AUC:",round(roctraingender1y$AUC,3)),
                 paste("Location AUC:",round(roctrainlocation1y$AUC,3)),
        
                 paste("T AUC:",round(roctraint1y$AUC,3)),
                 paste("Surgery AUC:",round(roctrainsurgery1y$AUC,3)),
                 paste("SRLNR AUC:",round(roctrainrlnr1y$AUC,3)),
                 paste("Radiation AUC:",round(roctrainradiation1y$AUC,3)),
                 paste("Chemotherapy AUC:",round(roctrainchemotherapy1y$AUC,3))),
       x.intersp=0.1, y.intersp=0.25,
       lty= 1 ,lwd= 2,col=c("red","green","blue","yellow","orange","purple","gray","pink","skyblue"),
       bty = "n",
       seg.len=0.5,cex=1)

# 训练集三年各变量比较图
plot(roctrain3y$FP,roctrain3y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctrainage3y$FP,roctrainage3y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctraingender3y$FP,roctraingender3y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
lines(roctrainlocation3y$FP,roctrainlocation3y$TP,type="l",col="yellow",xlim=c(0,1),ylim = c(0,1))
lines(roctraint3y$FP,roctraint3y$TP,type="l",col="orange",xlim=c(0,1),ylim = c(0,1))
lines(roctrainsurgery3y$FP,roctrainsurgery3y$TP,type="l",col="purple",xlim=c(0,1),ylim = c(0,1))
lines(roctrainrlnr3y$FP,roctrainrlnr3y$TP,type="l",col="gray",xlim=c(0,1),ylim = c(0,1))
lines(roctrainradiation3y$FP,roctrainradiation3y$TP,type="l",col="pink",xlim=c(0,1),ylim = c(0,1))
lines(roctrainchemotherapy3y$FP,roctrainchemotherapy3y$TP,type="l",col="skyblue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.6,c(paste("Nomogram AUC:",round(roctrain3y$AUC,3)),
                 paste("Age AUC:",round(roctrainage3y$AUC,3)),
                 paste("Gender AUC:",round(roctraingender3y$AUC,3)),
                 paste("Location AUC:",round(roctrainlocation3y$AUC,3)),

                 paste("T AUC:",round(roctraint3y$AUC,3)),
                 paste("Surgery AUC:",round(roctrainsurgery3y$AUC,3)),
                 paste("SRLNR AUC:",round(roctrainrlnr3y$AUC,3)),
                 paste("Radiation AUC:",round(roctrainradiation3y$AUC,3)),
                 paste("Chemotherapy AUC:",round(roctrainchemotherapy3y$AUC,3))),
       x.intersp=0.1, y.intersp=0.25,
       lty= 1 ,lwd= 2,col=c("red","green","blue","yellow","orange","purple","gray","pink","skyblue"),
       bty = "n",
       seg.len=0.5,cex=1)


# 训练集五年各变量比较图
plot(roctrain5y$FP,roctrain5y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctrainage5y$FP,roctrainage5y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctraingender5y$FP,roctraingender5y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
lines(roctrainlocation5y$FP,roctrainlocation5y$TP,type="l",col="yellow",xlim=c(0,1),ylim = c(0,1))

lines(roctraint5y$FP,roctraint5y$TP,type="l",col="orange",xlim=c(0,1),ylim = c(0,1))
lines(roctrainsurgery5y$FP,roctrainsurgery5y$TP,type="l",col="purple",xlim=c(0,1),ylim = c(0,1))
lines(roctrainrlnr5y$FP,roctrainrlnr5y$TP,type="l",col="gray",xlim=c(0,1),ylim = c(0,1))
lines(roctrainradiation5y$FP,roctrainradiation5y$TP,type="l",col="pink",xlim=c(0,1),ylim = c(0,1))
lines(roctrainchemotherapy5y$FP,roctrainchemotherapy5y$TP,type="l",col="skyblue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.6,c(paste("Nomogram AUC:",round(roctrain5y$AUC,3)),
                 paste("Age AUC:",round(roctrainage5y$AUC,3)),
                 paste("Gender AUC:",round(roctraingender5y$AUC,3)),
                 paste("Location AUC:",round(roctrainlocation5y$AUC,3)),
                 
                 paste("T AUC:",round(roctraint5y$AUC,3)),
                 paste("Surgery AUC:",round(roctrainsurgery5y$AUC,3)),
                 paste("SRLNR AUC:",round(roctrainrlnr5y$AUC,3)),
                 paste("Radiation AUC:",round(roctrainradiation5y$AUC,3)),
                 paste("Chemotherapy AUC:",round(roctrainchemotherapy5y$AUC,3))),
       x.intersp=0.1, y.intersp=0.25,
       lty= 1 ,lwd= 2,col=c("red","green","blue","yellow","orange","purple","gray","pink","skyblue"),
       bty = "n",
       seg.len=0.5,cex=1)


# 测试集各变量比较ROC
test$coxtrainagelp<-predict(coxtrainage,type="risk",newdata =test)
roctestage1y=survivalROC(Stime=test$Survival,
                         status=test$Status,
                         marker=test$coxtrainagelp,
                         predict.time=12,method="KM")
roctestage3y=survivalROC(Stime=test$Survival,
                         status=test$Status,
                         marker=test$coxtrainagelp,
                         predict.time=36,method="KM")
roctestage5y=survivalROC(Stime=test$Survival,
                         status=test$Status,
                         marker=test$coxtrainagelp,
                         predict.time=60,method="KM")


test$coxtraingenderlp<-predict(coxtraingender,type="risk",newdata =test)
roctestgender1y=survivalROC(Stime=test$Survival,
                            status=test$Status,
                            marker=test$coxtraingenderlp,
                            predict.time=12,method="KM")
roctestgender3y=survivalROC(Stime=test$Survival,
                            status=test$Status,
                            marker=test$coxtraingenderlp,
                            predict.time=36,method="KM")
roctestgender5y=survivalROC(Stime=test$Survival,
                            status=test$Status,
                            marker=test$coxtraingenderlp,
                            predict.time=60,method="KM")



test$coxtrainlocationlp<-predict(coxtrainlocation,type="risk",newdata =test)
roctestlocation1y=survivalROC(Stime=test$Survival,
                              status=test$Status,
                              marker=test$coxtrainlocationlp,
                              predict.time=12,method="KM")
roctestlocation3y=survivalROC(Stime=test$Survival,
                              status=test$Status,
                              marker=test$coxtrainlocationlp,
                              predict.time=36,method="KM")
roctestlocation5y=survivalROC(Stime=test$Survival,
                              status=test$Status,
                              marker=test$coxtrainlocationlp,
                              predict.time=60,method="KM")




test$coxtraintlp<-predict(coxtraint,type="risk",newdata =test)
roctestt1y=survivalROC(Stime=test$Survival,
                       status=test$Status,
                       marker=test$coxtraintlp,
                       predict.time=12,method="KM")
roctestt3y=survivalROC(Stime=test$Survival,
                       status=test$Status,
                       marker=test$coxtraintlp,
                       predict.time=36,method="KM")
roctestt5y=survivalROC(Stime=test$Survival,
                       status=test$Status,
                       marker=test$coxtraintlp,
                       predict.time=60,method="KM")



test$coxtrainsurgerylp<-predict(coxtrainsurgery,type="risk",newdata =test)
roctestsurgery1y=survivalROC(Stime=test$Survival,
                             status=test$Status,
                             marker=test$coxtrainsurgerylp,
                             predict.time=12,method="KM")
roctestsurgery3y=survivalROC(Stime=test$Survival,
                             status=test$Status,
                             marker=test$coxtrainsurgerylp,
                             predict.time=36,method="KM")
roctestsurgery5y=survivalROC(Stime=test$Survival,
                             status=test$Status,
                             marker=test$coxtrainsurgerylp,
                             predict.time=60,method="KM")



test$coxtrainrlnrlp<-predict(coxtrainrlnr,type="risk",newdata =test)
roctestrlnr1y=survivalROC(Stime=test$Survival,
                          status=test$Status,
                          marker=test$coxtrainrlnrlp,
                          predict.time=12,method="KM")
roctestrlnr3y=survivalROC(Stime=test$Survival,
                          status=test$Status,
                          marker=test$coxtrainrlnrlp,
                          predict.time=36,method="KM")
roctestrlnr5y=survivalROC(Stime=test$Survival,
                          status=test$Status,
                          marker=test$coxtrainrlnrlp,
                          predict.time=60,method="KM")



test$coxtrainradiationlp<-predict(coxtrainradiation,type="risk",newdata =test)
roctestradiation1y=survivalROC(Stime=test$Survival,
                               status=test$Status,
                               marker=test$coxtrainradiationlp,
                               predict.time=12,method="KM")
roctestradiation3y=survivalROC(Stime=test$Survival,
                               status=test$Status,
                               marker=test$coxtrainradiationlp,
                               predict.time=36,method="KM")
roctestradiation5y=survivalROC(Stime=test$Survival,
                               status=test$Status,
                               marker=test$coxtrainradiationlp,
                               predict.time=60,method="KM")



test$coxtrainchemotherapylp<-predict(coxtrainchemotherapy,type="risk",newdata =test)
roctestchemotherapy1y=survivalROC(Stime=test$Survival,
                                  status=test$Status,
                                  marker=test$coxtrainchemotherapylp,
                                  predict.time=12,method="KM")
roctestchemotherapy3y=survivalROC(Stime=test$Survival,
                                  status=test$Status,
                                  marker=test$coxtrainchemotherapylp,
                                  predict.time=36,method="KM")
roctestchemotherapy5y=survivalROC(Stime=test$Survival,
                                  status=test$Status,
                                  marker=test$coxtrainchemotherapylp,
                                  predict.time=60,method="KM")
# 测试集一年各变量比较图
plot(roctest1y$FP,roctest1y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctestage1y$FP,roctestage1y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctestgender1y$FP,roctestgender1y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
lines(roctestlocation1y$FP,roctestlocation1y$TP,type="l",col="yellow",xlim=c(0,1),ylim = c(0,1))

lines(roctestt1y$FP,roctestt1y$TP,type="l",col="orange",xlim=c(0,1),ylim = c(0,1))
lines(roctestsurgery1y$FP,roctestsurgery1y$TP,type="l",col="purple",xlim=c(0,1),ylim = c(0,1))
lines(roctestrlnr1y$FP,roctestrlnr1y$TP,type="l",col="gray",xlim=c(0,1),ylim = c(0,1))
lines(roctestradiation1y$FP,roctestradiation1y$TP,type="l",col="pink",xlim=c(0,1),ylim = c(0,1))
lines(roctestchemotherapy1y$FP,roctestchemotherapy1y$TP,type="l",col="skyblue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.6,c(paste("Nomogram AUC:",round(roctest1y$AUC,3)),
                 paste("Age AUC:",round(roctestage1y$AUC,3)),
                 paste("Gender AUC:",round(roctestgender1y$AUC,3)),
                 paste("Location AUC:",round(roctestlocation1y$AUC,3)),
                 
                 paste("T AUC:",round(roctestt1y$AUC,3)),
                 paste("Surgery AUC:",round(roctestsurgery1y$AUC,3)),
                 paste("SRLNR AUC:",round(roctestrlnr1y$AUC,3)),
                 paste("Radiation AUC:",round(roctestradiation1y$AUC,3)),
                 paste("Chemotherapy AUC:",round(roctestchemotherapy1y$AUC,3))),
       x.intersp=0.1, y.intersp=0.25,
       lty= 1 ,lwd= 2,col=c("red","green","blue","yellow","orange","purple","gray","pink","skyblue"),
       bty = "n",
       seg.len=0.5,cex=1)

# 测试集三年各变量比较图
plot(roctest3y$FP,roctest3y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctestage3y$FP,roctestage3y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctestgender3y$FP,roctestgender3y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
lines(roctestlocation3y$FP,roctestlocation3y$TP,type="l",col="yellow",xlim=c(0,1),ylim = c(0,1))

lines(roctestt3y$FP,roctestt3y$TP,type="l",col="orange",xlim=c(0,1),ylim = c(0,1))
lines(roctestsurgery3y$FP,roctestsurgery3y$TP,type="l",col="purple",xlim=c(0,1),ylim = c(0,1))
lines(roctestrlnr3y$FP,roctestrlnr3y$TP,type="l",col="gray",xlim=c(0,1),ylim = c(0,1))
lines(roctestradiation3y$FP,roctestradiation3y$TP,type="l",col="pink",xlim=c(0,1),ylim = c(0,1))
lines(roctestchemotherapy3y$FP,roctestchemotherapy3y$TP,type="l",col="skyblue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.6,c(paste("Nomogram AUC:",round(roctest3y$AUC,3)),
                 paste("Age AUC:",round(roctestage3y$AUC,3)),
                 paste("Gender AUC:",round(roctestgender3y$AUC,3)),
                 paste("Location AUC:",round(roctestlocation3y$AUC,3)),
             
                 paste("T AUC:",round(roctestt3y$AUC,3)),
                 paste("Surgery AUC:",round(roctestsurgery3y$AUC,3)),
                 paste("SRLNR AUC:",round(roctestrlnr3y$AUC,3)),
                 paste("Radiation AUC:",round(roctestradiation3y$AUC,3)),
                 paste("Chemotherapy AUC:",round(roctestchemotherapy3y$AUC,3))),
       x.intersp=0.1, y.intersp=0.25,
       lty= 1 ,lwd= 2,col=c("red","green","blue","yellow","orange","purple","gray","pink","skyblue"),
       bty = "n",
       seg.len=0.5,cex=1)


# 测试集五年各变量比较图
plot(roctest5y$FP,roctest5y$TP,type="l",col="red",xlim=c(0,1),ylim = c(0,1),
     xlab = paste("1-Specificity"),
     ylab="Sensitivity")
abline(0,1)
lines(roctestage5y$FP,roctestage5y$TP,type="l",col="green",xlim=c(0,1),ylim = c(0,1))
lines(roctestgender5y$FP,roctestgender5y$TP,type="l",col="blue",xlim=c(0,1),ylim = c(0,1))
lines(roctestlocation5y$FP,roctestlocation5y$TP,type="l",col="yellow",xlim=c(0,1),ylim = c(0,1))

lines(roctestt5y$FP,roctestt5y$TP,type="l",col="orange",xlim=c(0,1),ylim = c(0,1))
lines(roctestsurgery5y$FP,roctestsurgery5y$TP,type="l",col="purple",xlim=c(0,1),ylim = c(0,1))
lines(roctestrlnr5y$FP,roctestrlnr5y$TP,type="l",col="gray",xlim=c(0,1),ylim = c(0,1))
lines(roctestradiation5y$FP,roctestradiation5y$TP,type="l",col="pink",xlim=c(0,1),ylim = c(0,1))
lines(roctestchemotherapy5y$FP,roctestchemotherapy5y$TP,type="l",col="skyblue",xlim=c(0,1),ylim = c(0,1))
legend(0.5,0.6,c(paste("Nomogram AUC:",round(roctest5y$AUC,3)),
                 paste("Age AUC:",round(roctestage5y$AUC,3)),
                 paste("Gender AUC:",round(roctestgender5y$AUC,3)),
                 paste("Location AUC:",round(roctestlocation5y$AUC,3)),
                 
                 paste("T AUC:",round(roctestt5y$AUC,3)),
                 paste("Surgery AUC:",round(roctestsurgery5y$AUC,3)),
                 paste("SRLNR AUC:",round(roctestrlnr5y$AUC,3)),
                 paste("Radiation AUC:",round(roctestradiation5y$AUC,3)),
                 paste("Chemotherapy AUC:",round(roctestchemotherapy5y$AUC,3))),
       x.intersp=0.1, y.intersp=0.25,
       lty= 1 ,lwd= 2,col=c("red","green","blue","yellow","orange","purple","gray","pink","skyblue"),
       bty = "n",
       seg.len=0.5,cex=1)


cox_fittrain1y<-coxph(Surv(Survival,Status)~Age+Gender+Location+Grade+T+Surgery+SRLNR+Radiation+Chemotherapy,data=train,x=T,y=T)

cox_fits<-Score(list("1"=cox_fittrain1y),
                formula=Surv(Survival,Status) ~ 1,
                data=train,
                plots="calibration",
                conf.int=T,
                B=1000,
                M=50,
                times=12)
plotCalibration(cox_fits,
                xlab="Nomogram predicted probability of survival less than 12 months",
                ylab="Actual probability of survival less than 12 months",
               col = "blue",
               auc.in.legend = F,
               brier.in.legend = F)

cox_fits<-Score(list("1"=cox_fittrain1y),
                formula=Surv(Survival,Status) ~ 1,
                data=train,
                plots="calibration",
                conf.int=T,
                B=1000,
                M=50,
                times=36)
plotCalibration(cox_fits,
                xlab="Nomogram predicted probability of survival less than 36 months",
                ylab="Actual probability of survival less than 36 months",
                col = "blue",
                auc.in.legend = F,
                brier.in.legend = F)
cox_fits<-Score(list("1"=cox_fittrain1y),
                formula=Surv(Survival,Status) ~ 1,
                data=train,
                plots="calibration",
                conf.int=T,
                B=1000,
                M=50,
                times=60)
plotCalibration(cox_fits,
                xlab="Nomogram predicted probability of survival less than 60 months",
                ylab="Actual probability of survival less than 60 months",
                col = "blue",
                auc.in.legend = F,
                brier.in.legend = F)
cox_fits<-Score(list("1"=cox_fittrain1y),
                formula=Surv(Survival,Status) ~ 1,
                data=test,
                plots="calibration",
                conf.int=T,
                B=1000,
                M=50,
                times=12)
plotCalibration(cox_fits,
                xlab="Nomogram predicted probability of survival less than 12 months",
                ylab="Actual probability of survival less than 12 months",
                col = "blue",
                auc.in.legend = F,
                brier.in.legend = F)

cox_fits<-Score(list("1"=cox_fittrain1y),
                formula=Surv(Survival,Status) ~ 1,
                data=test,
                plots="calibration",
                conf.int=T,
                B=1000,
                M=50,
                times=36)
plotCalibration(cox_fits,
                xlab="Nomogram predicted probability of survival less than 36 months",
                ylab="Actual probability of survival less than 36 months",
                col = "blue",
                auc.in.legend = F,
                brier.in.legend = F)
cox_fits<-Score(list("1"=cox_fittrain1y),
                formula=Surv(Survival,Status) ~ 1,
                data=test,
                plots="calibration",
                conf.int=T,
                B=1000,
                M=50,
                times=60)
plotCalibration(cox_fits,
                xlab="Nomogram predicted probability of survival less than 60 months",
                ylab="Actual probability of survival less than 60 months",
                col = "blue",
                auc.in.legend = F,
                brier.in.legend = F)

train$risk<-ifelse(train$lp>median(train$lp),"high","low")

fit<-survfit(Surv(Survival,Status)~risk,data=train)
d<-data.frame(time=fit$time,
              n.risk=fit$n.risk,
              n.event=fit$n.event,
              surv=fit$surv,
              upper=fit$upper,
              lower=fit$lower)

ggsurvplot(fit,
           pval=TRUE,conf.int = TRUE,
           risk.table = TRUE,
           risk.table.col="strata",
           xlab="Time,months",
           linetype = "strata",
           surv.median.line = "hv",
           ggtheme=theme_bw(),
           palette = c("#E7B800","#2E9FDF"))

test$risk<-ifelse(test$lp>median(train$lp),"high","low")
fit1<-survfit(Surv(Survival,Status)~risk,data=test)
d1<-data.frame(time=fit1$time,
              n.risk=fit1$n.risk,
              n.event=fit1$n.event,
              surv=fit1$surv,
              upper=fit1$upper,
              lower=fit1$lower)
head(d)
ggsurvplot(fit1,
           pval=TRUE,conf.int = TRUE,
           risk.table = TRUE,
           risk.table.col="strata",
           xlab="Time,months",
           linetype = "strata",
           surv.median.line = "hv",
           ggtheme=theme_bw(),
           palette = c("#E7B800","#2E9FDF"))